Integration of Langevin Equations with Multiplicative Noise and the Viability of Field 

Theories for Absorbing Phase Transitions 
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EfRcient and accurate integration of stochastic (partial) differential equations with multiplicative 
noise can be obtained through a split-step scheme, which separates the integration of the deter- 
ministic part from that of the stochastic part, the latter being performed by sampling exactly the 
solution of the associated Fokker-Planck equation. We demonstrate the computational power of 
this method by applying it to most absorbing phase transitions for which Langevin equations have 
been proposed. This provides precise estimates of the associated scaling exponents, clarifying the 
classification of these nonequilibrium problems, and confirms or refutes some existing theories. 
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Stochastic differential equations are ubiquitous in the 
description of phenomena in the natural sciences and be- 
yond The coarse-graining of fast degrees of free- 
dom often leads to effective Langevin equations where 
the noise term involves the mesoscopic variables of in- 
terest in a multiplicative fashion. Examples range from 
nonlinear quantum optics, synchronization of oscillators, 
or wetting phenomena, to theoretical population dynam- 
ics studies and autocatalytic chemical reactions [ll |2, 13 ■ 

An important case in nonequilibrium statistical physics 
is the stochastic partial differential equation governing a 
single, positive, concentration field p = p{'r,t): 

dtpir, t)^D \j'^p + ap- bp^ + aV^?7(r, <), (1) 

where is a Gaussian (zero-mean) white noise (that is 
with correlations {ri{r,t)ri{r' ,t')) — 6{r — r')S{t ~ t')). 
For instance, for the reaction-diffusion process A 2A, 
2A 0, Eq. can be obtained in a variety of ways, ei- 
ther from phenomelogical considerations or through more 
rigorous transformations . Also named "Reggeon field 
theory" for historical reasons, Eq. describes the most 
prominent class of absorbing phase transitions (APT), 
the directed percolation (DP) class Indeed, inter- 
preted in the Ito (prepoint) sense, the unique, homoge- 
neous p — solution does not evolve: it is an absorbing 
state. Although a wealth of models have been found to 
exhibit a DP transition, this class does not encompass all 
possible cases, and the classification of APTs is currently 
a very active field 0, Q . Not only such an endeavor is 
of importance for conceptual reasons, but it should also 
yield a better understanding of the key ingredients which 
have impeded so far clear-cut experimental realizations 
of even the DP transition. 

Following this line of thought, stochastic equations 
similar to have been proposed as candidate field the- 
ories for related problems (see below). Their analyses 
are notoriously difficult, and mostly rely on the pertur- 
bative renormalization group machinery in the vicinity 



of the corresponding upper critical dimension, one of the 
few exceptions being a recent non-perturbative treatment 
of Eq. 1^ in 01 . Given this analytical bottleneck, it is 
tempting, with ever-improving numerical resources, to 
directly integrate such stochastic equations in order to 
check whether they at least exhibit the universal proper- 
ties they are supposed to represent. However, standard 
schemes either immediately run into severe difficulties. 
For instance, even for the zero-dimensional version of 
Eq. (Q, a first-order explicit Euler method, viz. p{t + 
At) p{t) + At[ap{t) - bp^ {t)] +ay/Atp{t)N{0, 1) , where 
A^(0, 1) is a normal random variate, will ineluctably pro- 
duce unphysical negative values for p{t + At), and all 
the more so when p —^ 0, the regime of interest for the 
APT. Another route, which would first trade the square- 
root noise for a less singular one through some change 
of variables (e.g., p ^ ip"^, or a Cole-Hopf transforma- 
tion p — > e^^"^), is also numerically unbearable since it 
generates pathological deterministic terms as the origi- 
nal variable p — > 0. Faced with this problem in the same 
context, Dickman proposed 0, somewhat ironically, to 
also discretize p, yielding a scheme consistent to the order 
0{^/Ai) in the limit At — > 0. This approach has been 
used with some success 0) S 13 ' but one can legitimately 
wonder to what extent one is truly simulating the origi- 
nal, continuous equation. In this respect, that the asso- 
ciated results are affected by the same long transients as 
those observed in microscopic models is also worrisome. 

In this Letter, elaborating upon a method pioneered 
by Pechenik and Levine in the somewhat distant context 
of front selection mechanisms in microscopic reaction- 
diffusion models we overcome the above hurdles. 
We first demonstrate the power of this approach on 
Eq. ^ before applying it to most related APTs for 
which a Langevin equation has been proposed, includ- 
ing the voter critical point with its two symmetric ab- 
sorbing states 0, Q|. Our results are particularly 
worthy in the context of the current debate about 
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APTs occurring when p is coupled to an auxiliary field 
■0: when tp is static and conserved (Manna sandpile 
model, conserved-DP, or fixed energy sandpiles class) 
Hi Qi Qi Q we obtain the best numerical esti- 

mates for the critical indices. When ip is conserved but 
diffuses El E US m m m 13, om results suggest 
that, at least in low spatial dimensions, the Langevin 
equations postulated or derived (approximately) as can- 
didate field theories are not viable. 

The idea underlying the approach of [lOl | (of which we 
become aware while this paper was upon completion) is 
a general and rather natural one, since it consists in in- 
tegrating the fast degrees of freedom. An Ito stochastic 
differential equation of the form ^ = f{p) + ag{p)ri{t) 
is dealt with the so-called operator-splitting scheme: the 
stochastic part ag(p)ri(t) is integrated first, not by using 
a Gaussian random number, but by directly sampling the 
time-dependent solution of the associated Fokker-Planck 
equation (FPE). Namely, one generates a random num- 
ber p* distributed according to the conditional transition 
probability density function (p.d.f.) Pioha..{p{t + At) = 
p*\p{t) — po}, and then use p* for evolving the deter- 
ministic part f(p) with any standard numerical method 
for ordinary differential equations. Since the integration 
of the stochastic part is accomplished through the ex- 
act solution of the FPE, which is first-order in time, the 
overall precision of the scheme, 0{At), is already signifi- 
cantly superior to that of a naive Euler method (anyhow 
flawed for Eq.lQJ), or to Dickman's approach. 

Now, for the square-root noise case, i.e., g{p) = the 
closed form solution P{p,t) = Proba.{p(t) = p\p{0) — 
Pq} of the associated FPE dtP{p,t) = ^d^p[pP{p,t)], 
has been known in the mathematical literature for more 
than half a century (see also 



a + Pp + <7y^r] using Eq.® with /? 



P{p,t)=S{p)e 



2(P0+P 




(2) 



(/i is the modified Bessel function of the first kind of 
order 1). When, further, the deterministic part is linear, 
i.e. f{p) ~ a + (3p, with a > 0, the exact conditional 
transition p.d.f. of the full equation ^ = a + Pp + a^i] 
has also been determined 1251: 
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2\^/^;;^) , (3) 



(/^ being a Bessel function of order p) where, to condense 



notations, we have set A 
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Y) , and p = —1 



The scheme we have used to integrate Eq.Q and its 
siblings relies on the latter results. After discretizing the 
Laplacian V^p over the 2d nearest-neighbors r + of 
site r on a d-dimensional hypercubic lattice of mesh size 
Ax, we first sample, between t and t + Ai, the solution 
of the FPE associated to each local linear equation ^ = 
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p(r + e„,i). 



(4) 



The value p* coming from the stochastic sampling step is, 
by construction, automatically non-negative, and serves 
as the initial condition for the remaining part of Eq.p|. 
i.e. dtp{v,t) ~ —bp'^{r,t), which can be trivially inte- 

= — Given that 6 > 0, 



grated to yield p{r, t + At) = Y+^^bSt ■ 
the non- negativity of p(r, t) will be preserved at all times 
if, initially, p{r) > everywhere, since a given by Eq.l^J 
will also be non-negative and Eq.® can be used. 

It remains to sample the above p.d.f., Eq.Q or Eq.©. 
Instead of using a table method, as the authors of 10], we 
remark that, with the help of the Taylor-series expansion 
of the Bessel function, Eq.® can be rewritten as 

' n! F n+u-l-l) ^ ' 

n—O ^ ^ 

In other words, one has the following mixture |2^ : 

p* = Gamma[^ + 1 + Poisson[A/9oe'^*]]/A, (6) 

where Prob.{Poisson[A/9oe'^*] = n} = i^lf^ — )_ f ^ 



r(^) 



This procedure 



and Prob.{Gamma[w] = v} = 

will reconstitute, on average, all the terms of Eq.© with 
their correct probability, and gives us — since standard 
and uniformly fast generators of Poisson and Gamma 
random numbers are available — a means of sampling 
in a "numerically exact" way these p.d.f. 0|. 

Typical results for Eq.Q in one dimension are shown 
in Fig. n along with data obtained using Dickman's 
method. Except for the (weak) linear stability require- 
ment coming from the discretized Laplacian, there is no 
limitation on At with the former method, so that the 
computational gain is of several orders of magnitude, to- 
gether with an unusually clean algebraic decay of {p), 
with an exponent 9 — 0.1595(2) matching to the fourth 
decimal the series-expansion result HQ. In fact, even if 
At — 0.25 for this run, the threshold a = adAt) is within 
one percent off its extrapolated limit value as At —^ 0, 
suggesting that the continuous limit of Eq. is already 
resolved. One of the reasons for the particular efficiency 
of this scheme even with such a large timestep is that it 
automatically takes into account, and in a self-adaptive 
fashion through the locally varying value of a (Eq.l^J), 
the strongly non-Gaussian modifications undergone by 
the instantaneous, conditional p.d.f., Eq.Q or Eq.©, as 
one gets closer and closer to the absorbing barrier. 

We now present some of our most salient results ob- 
tained for Langevin equations similar to deferring a 
more detailed account of our investigations to 28] . 

DP coupled to a non- conserved, non-diffusing field. To 
account for reaction-diffusion processes such as 2A — > SA, 
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FIG. 1: (Color online) Density decay (p) = {p{r,t))r ~ 
at criticality for Eq.Q in d = 1. Lower curve: using our 
scheme (At = 0.25, Ax = 1, a = Uc = 1.75623(2), b = 1, D ^ 
0.25, — 2, single run for a system of 2^^ sites with p = 1 
everywhere initially); a least-square fit gives 9 = 0.1595(2). 
Upper curve: using Dickman's method (similar conditions, 
but At = 10""^). Inset: Plateau of the local exponent. 



2A — > A where single particles do not move (such as 
the prototypical pair-contact process [i^) and which 
thus possess infinitely-many absorbing states, it has 
been proposed that Eq.(^) be supplemented by the 
non-Markovian term cexp[— w J^dsp{r,s)] to account 
for the memory effect introduced by immobile parti- 
cles [sOl- The impact of this term is however unclear, 
with early simulations |^ using Dickman's method sug- 
gesting that continuously-varying spreading exponents 
arise, in agreement with results obtained on microscopic 
models [29j, but in contradiction with the study of 
infinite- memory spreading processes jSlJ, which support 
stretched-exponential behavior. Simulations with our 
scheme, in one dimension, reveal power-laws for small |c|, 
but curvature appears at late times for large, negative c 
values (Fig.l^L). To be fully conclusive, these results will 
have to be improved by using enrichment methods en- 
abling to explore rare events, but they already indicate 
that the conclusions of ^] probably hold asymptotically. 
We finally mention that in two dimensions (and for c > 0, 
6 = 0) we obtain dynamical percolation spreading expo- 
nents as predicted by the standard theory [soj . 

DP coupled to a conserved, non- diffusing field. Reac- 
tion processes such as A + B 2A, A ^ B where A 
particles diffuse and B are static, are similar to the case 
above, but the number of particles is (locally) conserved 
[l^ . This conservation law leads to couple p to a con- 
served field (j) in the following system (1^ J^] : 

dtp = DV'^p + ap-bp^ +ujp(j) + ij^ri{v,t) , , 



Microscopic models leading to (O also include so-called 
fixed energy sandpiles such as the Manna model, estab- 
lishing a link between APTs and self-organized critical- 
ity [l6j. The conservation law influences even the static 
exponents but definite estimates are currently not avail- 
able (see and references therein). Data from micro- 
scopic models, as well as simulations of ^ using Dick- 
man's method are plagued by long transients/corrections 
to scaling. Our scheme leads, again, to clean power-laws 
which provide us with the best estimates for the scal- 
ing exponents of this class of APT. In Fig. we show 
a typical result for critical decay m d = 1, leading to 
6 = 0.125(2), unambiguously distinct from the DP value 
0.1595(1). Critical decay exponents obtained in higher 
dimensions, 6*2^ = 0.509(5) and 03d = 0.81(1) 22], differ 
also significantly from their DP counterparts. 
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FIG. 2: (Color online) (a): Survival probability of an ini- 
tial seed at r = evolving under Eq.Q with the extra 
term cexp[—w f^^ dsp(r, s)]. Same parameters as in Fig. 
and various c values (about 10^ trials), (b): Same as Fig. Q 
but for Eqs.JTJ. Lower curve: with our scheme {At — 0.1, 
Ax = 1, a = ttc = 0.86455(5), 6 = lj = 1, D = Dp= 0.25, 
— 2, p = <j> = 1 everywhere initially); A least-square fit 
gives 6 = 0.124(1). Upper curve: with Dickman's method 
{At = 0.0025). 

DP coupled to a conserved, diffusing field. If, for the 
reaction processes above, both species are diffusing, the 
situation changes again, if only because one has now 
a single, dynamic absorbing state (where B particles 
diffuse in the absence of As) . This case was studied 
both analytically [llllUla HI H^l and numerically 
[2^ 01 with continuous APT predicted and observed for 
< < p , but with conflicting estimates of scaling ex- 
ponents "[i^, IS- The corresponding Langevin equation 
is usually cast |20l EH as Eqs.Q complemented by the 
self-diffusion of the auxiliary fleld and a conserved noise 
term. Performing with our scheme critical decay experi- 
ments in spatial dimensions d — 1,2 we flnd the exponent 
6 to be undistinguishable from the DP values. Because 
this differs from both analytical predictions [l^ 
and estimates from microscopic models this indi- 
cates that the truncation of the full action of the fleld 
theory needed to arrive at the corresponding Langevin 
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equations is not legitimate. 
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FIG. 3: (Color online) Voter coarsening process for a contin- 
uous field evolving under Eq.(|SJ (Ax = 1, At = 0.25, D = 1, 
— 0.5, zero-mean random initial conditions) Left: snap- 
shot of p at t = 2500 (512^ sites). Right: decay of interface 
density (4096^ sites) pi = \ - (p(r, t)p(r + e„, t)) cx 1/lnt. 



The voter critical point. The universality class of the 
voter model is characterized by two symmetric absorbing 
states ^3 . The following field theory has been proposed 
— but never tested — to describe its critical point ,lL.2l| : 



dtp ^DW^p + o-vT~~^?7(r, t) 



(8) 



The FPE associated to the sole stochastic part can be 
solved through an eigenfunction expansion, leading to 
a complicated expression for the conditional transition 
p.d.f., involving a continuous part and two delta peaks at 
the barriers p = ±1 0,|23- Although this distribution 
can be sampled '2^, it is both much simpler and more 
efficient to replace in the noise term the piece yl — by 
Q{p)\/1 — p+ (p ^ ^p), thereby taking into account just 
the closest DP barrier. This way our scheme can be ap- 
plied and, one observes in two dimensions phase ordering 
patterns typical of the marginal voter coarsening process, 
and, for the first time with continuous variables, the ex- 
pected 1/ \nt decay of the density of interfaces (Fig.l^J. 

This completes our (not-exhaustive) inspection of 
Langevin equations proposed as field theories of absorb- 
ing phase transitions. Pending more comprehensive stud- 
ies (higher dimensions, other scaling exponents), the re- 
sults already obtained demonstrate that the method pre- 
sented above enables faithful and efficient simulations of 
such stochastic equations. This approach will remain 
particularly useful as long as no major analytical progress 
is made, and also to test future theoretical predictions. 
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